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Abstract 

We consider reflector imaging in a weakly random waveguide. We address the sit- 
uation in which the source is farther from the reflector to be imaged than the energy 
cquipartition distance, but the receiver array is closer to the reflector to be imaged than 
the energy equipartition distance. As a consequence, the reflector is illuminated by a 
partially coherent field and the signals recorded by the receiver array are noisy. This 
paper shows that migration of the recorded signals cannot give a good image, but an 
appropriate migration of the cross correlations of the recorded signals can give a very 
good image. The resolution and stability analysis of this original functional shows that 
the reflector can be localized with an accuracy of the order of the wavelength even 
when the receiver array has small aperture, and that broadband sources are necessary 
to ensure statistical stability, whatever the aperture of the array. 

1 Introduction 

Sensor array imaging in a scattering medium is limited because coherent signals recorded at 
the source-receiver array and coming from a reflector to be imaged are dominated by inco- 
herent signals coming from multiple scattering by the medium. For instance, in a randomly 
perturbed waveguide, it is known that the field becomes completely incoherent when the 
propagation distance becomes larger than the equipartition distance, which corresponds to 
the distance beyond which the source energy has been shared equally among all the propa- 
gating modes |10t Chapter 20]. As we will see, if the distance between the source-receiver 
array and the reflector is larger than the equipartition distance, then classical migration of 
the signals recorded at the array cannot give a good image. 

Sources can be expensive or difficult to implement but receivers can be cheap and easy 
to implement, so an imaging problem in which there are a few sources (all of them being 
far from the reflector) and many receivers (some of them being close to the reflector) is 
of theoretical and practical interest. If there is a unique source far from the reflector 
(farther than the equipartition distance) and if the receiver array is close to the reflector 
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(closer than the equipartition distance), then classical migration of the recorded signals 
fails again. This was shown is various contexts and we will show it again in the waveguide 
geometry. However, in such a situation, another kind of migration can be used: from the 
work devoted to coherent interferometry imaging [U [5l [6j 18] and ambient noise imaging 
[9|ll2 1 [T6 l ll7 1 ll8j. it is known that migration of cross correlations of noisy signals can be more 
stable than migration of the signals themselves. The migration of cross correlations of noisy 
signals recorded by auxiliary passive arrays was proposed by [3] in geophysical contexts and 
analyzed recently in randomly scattering open media in |14j . and we would like to address 
the same problem in the waveguide geometry. Indeed the number of propagating modes is 
finite in the waveguide geometry so that the statistical behavior of partially coherent fields 
in random waveguides is very different from the open medium case I10j . In our paper, 
we show that, if a receiver array can be placed close to the reflector to be imaged, then the 
cross correlations of the incoherent signals on this array can be used to image the reflector. 
We will give a detailed resolution and stability analysis. We will show that the statistical 
stability requires a broadband source and that good resolution and stability properties do 
not require the receiver array to span the whole cross section of the waveguide, which is an 
effect specific to the waveguide geometry. 

The paper is organized as follows. In section O we review the mathematical background 
of the imaging problem in a random waveguide. In section [3] we describe and analyze the 
classical migration functional using the recorded signals and show that it cannot give a 
good image when the propagation distance is beyond the energy equipartition distance. In 
section [J we introduce the correlation-based imaging functional; it has two versions which 
correspond to the time-harmonic case and the broadband case. In section [5j we analyze 
the resolution of the proposed imaging functionals. Detailed analyses are provided for full 
aperture and limited aperture arrays. These results are based on the statistical average 
of the imaging functionals. The variances of these functionals are very important as well 
because they determine the statistical stability of the imaging functionals. In section [H we 
study the variances of the imaging functionals. Some concluding remarks are listed at the 
end of the paper. 

2 Mathematical Formulation of the Imaging Problem 

2.1 The ideal waveguide 

We consider linear scalar (acoustic) waves propagating in a two-dimensional space. The 
governing equation is 



Here p is the scalar field (acoustic pressure); cq is the speed of propagation in the medium 
(sound speed); F(t,x) models the forcing term. We consider a waveguide geometry, and we 
decompose the spatial variable x as (x, z). That is, z 6 1 is along the axis of the waveguide 
while x G T> denotes the transverse coordinate, and T> = (0, a) is the transverse section of 
the waveguide. We assume that the forcing term is localized in the plane z = 0: 






(2.2) 
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where x sourcc = (x s ,0) for some x s E V. We assume that the medium is quiescent before 
the pulse emission, that is 

p(i,x)=0, i<0. (2.3) 
We consider Dirichlet boundary conditions a the boundary of the waveguide: 

p(t,x) = 0, xeaPxl = {0,a}xi (2.4) 

Using the Fourier method, the scalar field can be written as a superposition of waveguide 
modes. A waveguide mode is a time-harmonic wave of the form p(uj,x.)e~ lLUt with frequency 
cj, where p satisfies the time-harmonic form of the wave equation (12. ip without a source 
term: 

d 2 p(oj, x, z) + A±p(oj, x, z) + k 2 (uj)p(uj, x, z) = 0. (2-5) 

Here, Aj_ = d 2 is the transverse Laplace operator in the transverse section T> with Dirich- 
let boundary conditions; k{uj) = uj/cq is the homogeneous wavenumber. Consequently, 
(|2.5p can be solved using the eigenmodes of Aj_, that is, using the orthonormal basis 
{'>.}{■>■)). i 1.2,- of L 2 (D) given by 

-A ± 4> j (x) = Xjcpjix), x£V. (2.6) 

The eigenvalues are simple, satisfying < Ai < A2 < • • • . The eigenvalues and eigenvectors 
are given by 

\/2 7rjx tt^J - ^ 
<t>j( x ) = "^ sm ( )> X j = — 2~- ( 2 - 7 ) 

Using the method of separating variables on (|2.5|) . we see that the waveguide mode p(oj, x, z) 
can be further written as superposition of pj(oo,x,z) = 4>j(x)e :izl ^ : >^ z , where 

f3](u}) = k 2 (ta)-\ j . (2.8) 

For a given frequency cj, there exists a unique integer N(uj) such that Xn^) — k 2 (to) < 

N(u>)= — . (2.9) 

L7TCoJ 

Here and in the sequel, |^J means the integer part of a real number b. The modes 
{pj(uj, x, z) = 4>j( x ) e±ll3j( ' UJ ^ z }j=i,...,N(ui) are propagating waveguide modes and {/^ (w)}./^...^^) 
are called the modal wavenumbers. On the other hand, {pj(uj, x, z) = 4'j(x)e ± ^^ ul ^ z }j >N ^ 
are evanescent modes because they decay as z goes to =Foo. 

2.2 The randomly perturbed waveguide 

From now on we assume that the waveguide is randomly perturbed and the scalar field 
satisfies the perturbed wave equation 

1 2 

Ap(i,x)-— ^|(i,x) =F(t,x), (2.10) 

where c(x) is the randomly heterogeneous speed of propagation of the medium. We consider 
the case where the typical amplitude of the fluctuations of the speed of propagation is small, 
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which we call the weakly random regime. When the correlation length of the fluctuations 
is of the same order as the typical wavelength the interactions between the waves and the 
random medium become nontrivial. Due to the small amplitude of the fluctuations, however, 
the effect of the interaction becomes important only after a long propagating distance. 

More exactly we assume that a randomly heterogeneous section in z S [0, Lq] is sand- 
wiched in between two homogeneous waveguides: The speed of propagation is of the form 

1 f X(l + eu(x,z)) if (x, z) £ [0, a] x [0, Lq], 

-o (2-11) 



c 2 (x,z) I otherwise 

Here, v is a mean-zero, stationary and ergodic random processes with respect to the axis 
coordinate z. It is assumed to satisfy strong mixing conditions in z. The relative amplitude 
of the fluctuations of the speed of propagation is denoted by e. We assume that the cor- 
relation length of the random perturbation is of the same order as the typical wavelength 
A(wo) = 2itcq/ujq = 2ir/k{(jJo), for ujq the central frequency of the source. We assume that 
the propagation distance Lq is much larger than the typical wavelength. We will see that 
the interesting regime is when the ratio X/Lq is of order e 2 , so we introduce the normalized 
propagation distance Lq: 

j _Lq 

In this regime the cumulative interaction of the scalar wave with the small fluctuations of 
the speed of propagation becomes of order one. 

For a fixed frequency cj, the Fourier transformed scalar field p(oj,x,z) defined by 

iu)t , 



p(oj,x,z) = J p(t,x, z)e luJ dt 
satisfies the equation 

d 2 z p(u), x, z) + A±p(ui, x, z) + k 2 (u)[l + ev{x, z)]p(u>, x, z) = f(oj)5(x — x s )5(z). (2.12) 

To solve this equation, we make the following two simplifications that are justified in |10t 
Chapter 20] or p]. 

Ignoring the evanescent modes. First, we only consider the propagating modes: 

N(u) 

p(u,x,z) = <t>j{x)Pj(v,z). (2.13) 
j'=i 

This is valid because we are mainly concerned with the scalar field for z 3> 1 and the 
evanescent modes decay exponentially fast. Furthermore, we parameterize the complex 
mode amplitude pj(uj,z) by the amplitudes of its right- and left-going components. Let 
dj(oj, z) and bj(u, z) be the amplitudes of these components, defined by 



pj = -= (a je ^ + bje-^) , = ( - b 3 e~^ ) . 
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Using these representations, one obtains a system of ordinary differential equations (ODEs) 
for {aj,bj} [10], Section 20.2.4] or [TTJ Section 3.1]. The coefficients of the system depends 
on the integrated quantities of the form 



Cji(z) = / <j)j(x)<f)i(x)v(x,z)dx. (2.15) 
Jv 

This system of ODEs is closed by the boundary conditions at z = where the source F is 
imposed and at z = Lq = Lq/s 2 where there is no left-going component. 

Forward scattering approximation. Second, we neglect the left-going (backward) propa- 
gating mode, assuming that they do not interact with the right-going ones. This is valid in 
the limit £ — > when the second-order moments of v satisfy certain conditions |10l Section 
20.2.6] (or [TTJ Section 3.3]). In this case, the rescaled amplitude 

a){z) = aj(uj,z/e 2 ) 

of the right-going wave satisfies 

£ - (?) ^ <"«> 

where a e denotes the iV(u;)-dimensional vector (af , . . . , aL^,)' and is a N(uS) x N(uS) 
complex matrix with components 

( o) _ ik 2 Cji(z) m - 0i)z , , 

Define the propagator matrix T £ (cj, z) to be the fundamental solution of the system (|2.16p . 
i.e., 

dT e 1 / z \ 

-(,,z) = -HWyT E (,,4 (2.18) 

with T £ (u;,0) = I. Then a e (uj,z) = T £ (u,z)a £ (uj,0), where the initial amplitude a £ (uj,0) 
is determined by the source F. In fact, integrating (I2.12p across he plane z = and using 
(I2.12I2.14|) . we find that 

ai(u,0) = L^ f(u)<f> t {x B ), 1 = 1,..., N(u>). 

2ly//3l(u) 

Consider an array of receivers located in the plane z = L of the random waveguide 
section, where 

L = 4, (2-19) 

and < L < Lq. Let T%(uj) be the j7-entry of the propagator matrix T £ (cj,L) at L. 
Then T% is the rate of the initial /-mode being scattered to the j-mode in the plane with 
z = L = L/e 2 of the random waveguide section. In particular, we have 

JV(w) JV(w) 

aj(u,L) = a £ (co,L) = £ T^af^O) = ^ — 7 _T^( W )/(a;)^(x s ). (2.20) 

1=1 i=i 2«V«(w) 
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Figure 1: Schematic of the imaging problem. A point source (circle) in the plane z = 
emits a short pulse that propagates through the random waveguide. The target (square) in 
the plane z = z r is a reflector. The receiver array (triangles) in the plane z = L records the 
signals. 



Since the typical amplitude of the random fluctuations of the speed of propagation is of 
order e, the amplitudes of the right-going waves do not change for z > L: 

aj{u,z) = CLj(u, L) + 0(s) 

provided z — L = O(l). Using these expressions, (I2.13P and (12.141) . we obtain the expression 
of the scalar field at z > L with z — L = 0(1): 

p(u,x,z) = £ 2j { ( w )/( w )^(x)^(a:8)e^*. (2.21) 



2.3 Modeling the point reflector 

In the imaging problem to be investigated (see Figure [I]), the goal is to locate a point 
reflector centered at x r = (x T ,z r ) from signals recorded at the receiver array in the plane 
z = L = L/e 2 . The reflector is supposed to be at a relatively small distance (compared to 
e -2 ), that is to say 

z r = ^ + (2-22) 

The reflector can be modeled as a local change in the density and/or the bulk modulus of 
the medium, so that the sound speed is locally modified as 

sfkrg + g 1 *^- <2 - 23) 

where f2 r is a small domain around x r := (x r , z r ) which represents the center of the reflector; 
c r is a parameter characterizing the contrast of the reflector. With this modification, the 
right-hand side of (|2.5p should have an additional term — (w/c r ) 2 ln r p(a;, x, z). We assume 
that the diameter of the scattering region f2 r is small compared to the typical wavelength 
and that the velocity contrast is such that cr r := c~ 2 |r2 r | satisfies a r <C 1. Then we can 
model the scattering region by a point reflector 

-ylQ t (x,z) « 0- r 5(x - X r ). 
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Born approximation. The above setting allows us to solve the scalar field with the 
presence of the point reflector using the Born approximation for the reflector. Given a fixed 
frequency lj, we have 

p(u, x, z) « p p (uj, x, z) + p s (u, x, z). (2.24) 

Here, p p is the primary field induced by the source F propagating through the random 
waveguide and computed in the previous section (Eq. (|2.2ip ). and p s is the secondary field, 
that is the first-order scattered field due to the additional source — w 2 <7 r 5(x — x r )j5 p at the 
reflector: 

Ap s (u), x, z) + k 2 (u) [l + eu(x, z)]p s (uj, x, z) = — aj 2 <7 r <5(x — x r )j5 p (u;, x T , z r ). (2.25) 

Note that in the Born approximation one replaces the full wave field at the reflector by the 
primary field. 

The primary field is solved exactly as in the previous section. Summarizing the results 
there, one obtains that for z > L and z — L = 0(1), 

P p (oj,x,z) = Y, J^> T^u)^(x)M^)e^ z . (2.26) 

The secondary field satisfies (|2.25j) . Again, we solve this equation using the orthonormal 

basis {^j(i)}j = i m w ) and we ignore the evanescent modes. Since the reflector is within 

a distance of order one from the receiver array, for L < z < z T , the propagator matrix 
can be approximated by the identity matrix up to a term of order e, and we only need to 
decompose the secondary source at the reflector into waveguide modes. Using (|2.25p . the 
decomposition (|2.14p and the fact that there is no left-going wave from z > z T , we find that 

. 2 

KjiU'Z) = n /^ fe^rfe^^rJr). (2-27) 



We note that there is no right-going secondary wave because we do not consider back- 
scattering of the left-going secondary wave near the receivers. Finally, recall the expression 
of the primary field at the reflector (|2.26|) . we obtain for z £ [L, z T ) that 



N(u) 

p s (u,x,z) = ]T^.(x)^^>,z)e-^-^ 

Pi 



j=l V ^ 

3,1^=1 WjVPmVPl 

Since c r <C 1 is a very small number, the secondary field field is dominated by that due to 
the primary field. 



3 Migration-Based Imaging Functional 

In this section, we introduce the classical imaging functional to localize the point reflector 
using the scalar (pressure) field recorded at the receiver array at z = L. This imaging 
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functional is based on the migration of the array data to a search point (x^S 8 ). Our goal 
is to show that classical Kirchhoff migration functional does not give a good image when 
the medium between the source at z = and the receiver array at z = L is scattering. 
The data of scalar (pressure) field recorded by the receivers are 

{p(t,x,L) | t G R,x G V). 

Note that we consider in this section the full aperture case: the receivers span the whole 
cross section of the waveguide and they record data at all time. We consider the frequency- 
and mode-dependent data 



pj{ui,z = L) = J J p{t,x,z = L)(f>j(x)dxe tujt dt. 

According to the analysis carried out Section 12.31 it can be decomposed as 

Pj(u),z = L) =p pj (u,z = L) +p sj (u,z = L). 
Prom (12.26|) and (|2.28j) the primary and secondary contributions are 

p pj (uj,z = L) = /7r e ? /b^Wb), (3.1) 

p sj {uj,z = L) = — -^- < /) i (x r )e 4/3 ^ r g(w,2; r ,2 r ), (3.2) 

Pj[UJ) 



with 



N 



which can be interpreted as an illumination of the reflector. The secondary contribution 
p s j contains the information about the reflector, and its form (|3.2p motivates the one of the 
Kirchhoff migration imaging functional: 

N(lu) 

Zkm(x s ,z s ) = - — £ P j (u)cl> j (x s )e- i ^ zS p j (uj,z = L)du>, (3.3) 

where the search point is (x s , 2 s ) with 2 s = L + z s . 

A simple case is when the source term is time-harmonic, i.e. F(t, x) = 5(x — x source )/(i) 
with f(t) = e~ iuJot and 

f(oj) = 2it5(uj - w ). 

Then the data set is reduced to {pj(ujq, z = L), j = 1, . . . , N(uo)} and the Kirchhoff migra- 
tion functional has the form 

JV(w ) 

1km(x S ,z s ) = -— V P j (u; )ci )j (x S )e- i ^> S p j (LO ,z = L) (3.4) 

We need to compute the mean of the imaging functional in order to characterize its res- 
olution properties and its variance in order to characterize its stability properties. These 
statistical moments depend on the moments of the propagator matrix which were studied 
in |1CH Propositions 20.6 and 20.8] or [11, Propositions 6.1 and 6.3]. 
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Proposition 3.1. The first-order moments of the transmission coefficients have limits as 
e — > 0, which are given by 

E[Tj^)]^0, ifj^l, (3.5) 

E[2£(w)] ^> e - D ^) L , otherwise. (3.6) 

T/ie second-order moments of the transmission coefficients have limits as e — > 0, which are 
given by 



E[T£( W )l|(u,)] ^> e-VM^ i/j ^ J, (3.7) 

E[^M75H]^>7- (0 (a;), (3.8) 

E[2j z (o;)7^(£j)] 0, otherwise. (3.9) 

T/ie functions Tj 1 \uj,L) are the solutions of the system of linear equations 



d z ^,-jnV~>\'n '3 y '3 

™¥=3 



V rgn (r« - 7f ) , 7?V, * = o) = s sl . (3.10) 



(c) 

The positive coefficients Dj and Qji and the matrix r^ n depend on the correlation function 
of the random process v. Furthermore, we have 



sup|EpT? { (w)]| < Ce~ L/L ^, sup 
3,1 3,1 



< Ce~ L/Lc ^, (3.11) 



where L cqu ip is the equipartition distance for the mean mode powers introduced at the end 
of Section 20.3.3 in (or at the end of Section 4.2 in f77]/). 

The results on the first-order moments describe how the wave loses its coherence as it 
propagates in the random waveguide. The results on the second-order moments describe 
how the wave energy becomes equipartitioned on the waveguide modes. 

When L is larger than the energy equipartition length -L e quipj then the first-order mo- 
ments of the transmission coefficients are vanishing. Based on this observation, we have 

E[Z K m(* S ^ S )] ~0. (3.12) 

It turns out that the fluctuations of the imaging functional are much larger than its mean. 
This can be seen by studying the standard deviation of the imaging functional. When L is 
larger than the energy equipartition length L e q U ip, then the second-order moments of the 
transmission coefficients are vanishing except E^Tj/j 2 ] which converge to 1/N. Based on 
this observation, the second-order moment of the imaging functional for a time-harmonic 
source is: 



Z , Xr , Z r 



E[|X KM (x s ,z s )| 2 ] = |±g l(a S) + (^1) 

+ (^J^) 9m(^(x S , -z s ; x r , -z v )V(x S , z S ; x r , -z T ))] , (3.13) 
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where, for any integer j, we have defined 



1 



N 




(3.14) 



n=l 



1 



N 



^0 n (x r )^ n (x s ) e ^»( 



(3.15) 



n=l 



The first term in the right-hand side of (|3. 13|) is the contribution of the primary field. 
The second term is the contribution of the secondary field. The third term is a crossed 
contribution. 

These results show that, when the waveguide is randomly perturbed and long enough 
(longer than the equipartition distance), then the illumination of the reflector becomes 
incoherent and Kirchhoff migration, which is based on coherent effects, gives a completely 
unstable and noisy image. 

The analysis is complete in the time-harmonic case. The analysis of the broadband case 
(when the support of the source spectrum is not reduced to a single carrier frequency) goes 
along the same line although it is necessary to use the asymptotic expressions of the two- 
frequency second-order moments of the transmission coefficients (see [10, Proposition 20.7] 
or jlH Proposition 6.3]): due to the loss of coherence, the mean of the imaging functional 
is zero while its variance is not. 

4 Correlation-Based Imaging Functionals 

In this section, we introduce a new imaging functional to localize the point reflector using 
the scalar field recorded at the receiver array at z = L. This functional is based on the 
correlation functions of the recorded signals, which we stress in the first subsection. 

4.1 Correlation of the scalar field 

Let A denote the positions of the receivers in the plane z = L. The data of scalar (pressure) 
field recorded by the receivers are {p e (t,x,L) \ t 6 R, x £ .4,}. For simplicity, we have 
assumed that the receivers record data at all time. From these data one can form the cross 
correlation of the recorded field: 



Using the decomposition p = p p + p s in (|2.24|) . we can decompose the above cross cor- 
relation function into four parts. Let C pp denote the cross correlation between the primary 




(4.1) 



In Fourier domain, it has the form: 




(4.2) 
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fields at the two receivers. Thanks to the formula (I2.26p . it admits the expression 
1 f N(UJ) 1 _ 

Cpp(r,X 1 ,X 2 )=— Yl loo o o 7 ^ T jl( U; ) T mn(^)\f(^)\ 2c Pj( x l) ( Pm(x2) 

87T J j,l,m^=l VPlPnPjPmiu) (4.3) 

Let C ps denote the cross correlation between the primary field at the first receiver with 
the secondary field at the second receiver. Recall that the secondary field contains infor- 
mation about the waves emitted from the reflector at x r = (x r ,z r ), with z r = L + z r . Due 
to (I2.26P and (I2.28p . it admits the expression 

f NH iu 2 a 
C ps (t,x 1 ,x 2 ) = 77 o 7Wn a ? , T jl^)T^ n {u)\f{u)\ 2 (t>j{xi)<l> q {x2) 

J qj^n^^^VfrPnPjPrnH (4.4) 

/ (x s )0n(xs)^(x r )<A m (x r )e 4 ^-^) Z e i ^ + ^)^e-^^. 

Similarly, let C sp denote the cross correlation between the secondary field at the first 
receiver with the primary field at the second receiver. One verifies that 

C sp (r, Xl ,x 2 ) = / £ - ^ ^ ==JJ f H7^ n (a;)|/(a;)| 2 ^(gi)0m(x2) /a . 

J q,j,l7?,n=l ™nMPlPnPjPm\u) (4-5) 

^Kx s )^ n (x s )^(x r )^(x r ) e i ^-^) I e- i ^ + ^> r e-^da;. 

Finally, the cross correlation between the secondary fields at the two receivers is much 
smaller than those above and its contribution is ignored. We neglected also the contributions 
from the error terms of the decomposition (|2.24p . These are justified because cr r <§C 1 
consistently with the Born approximation. 

Recall that the source in the acoustic model is due to the force F(t, x) = f(t)5(x— x sourcc ) 
where x source = (x s , 0) indicates the location of the source. In the rest of the paper, we will 
consider two special cases as follows. 

4.1.1 Cross correlation for broadband pulse 

We first consider the case where the source is given by F(t,x) = f(t)5(x — x source ) with 

f(t) = He a t)e~^ 1 . (4.6) 
Here, loq is the carrier frequency. In the Fourier domain, we have 

Here, /o is assumed to be a function with compact support or fast decay. 

When a > 2 the bandwidth has no effect and the situation is equivalent to the time- 
harmonic case that we address in the next section. 
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When a = (0, 2), the pulse is said to be broadband and the bandwidth plays a role in the 
propagation in the waveguide for a propagation distance of the order of L = L/e 2 . Although 
the analysis can be carried out in general, we restrict ourselves to the case a € (1, 2) because 
when a < 1, the number of propagating modes N(uj) varies with oj over the bandwidth and 
the analysis is a little bit more delicate. Nevertheless, the overall picture does not change 
in the latter case. 

Henceforth, a is a fixed number in the interval (1,2). Let u = luq + e a h. Then 
/>) = ^fo(h), 2J { (a;) = T^uo + e a h), 
in terms of the new variable h. Further, we have the following Taylor expansions 
Mu,) = fJ j +e»# i h + o(e») i / = 1 +0(0- 

yJPjP m PlPn\u) ^PjPmPlPn 

Here, is the derivative of fij at the carrier frequency ojq; further, the reduced wavenumber 
j3j is also evaluated at ojq. Using these formulas, the cross correlation functions become 

j,l,m,n=l V PjyrnPlh'n 

^(x s ) < A„(x s )e i [ /3m ^ 0+£Q/l )-ft^ 0+£Q ' l )] Z e- i(aJo+£Q/l)r d/ i , 



<M*s)<M^ m (*rK [ ^ 0+ ^ 

C *v~l E T^fe wPw wWi^ + ea/i ) T -(^o + 5 Q / i )^(x 1 )0 m (x 2 )^(x s ) 

3 q,j,l,m,n=l W7TP 1 £ VPjPmPlPn 



4.1.2 Cross correlation with time-harmonic source 

A simple case is when the source term is time-harmonic, i.e. F(t,x) = 5(x — x sourcc )/(i) 
with f(t) = e~ iuJot and 

f(oj) = 2tt5(uj - u ). 

Formally, one substitutes the above relation into the expressions (|4.3|) . (|4.4p and (|4.5|) 
and observes that the cross correlations are simplified to 

Cp P (r,xi,x 2 ) =- V ■===== Tf l (u )7^ in (u )4>j(xi)^ m (x2) 
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0« (z s ) <t> n {x s ) </> g (x r ) 4>j (x r ) e l 
More precisely, when the source is time-harmonic, the solution is of the form p(t, x, z) 



s i(J} m -Pj)L e -i(f) q +Pj)zi e -iuoT 



p(u)Q,x, z)e The definition of the correlation function should be modified to 

1 fT „ „ . _ „ 

C(t, Xl ,x 2 ) := - / p(t, xi,L)p(t + r, x 2 , L)dt = e~ iUJ0T p(xi,L)p(x 2 , L). (4.7) 
1 Jo 

Using the decomposition (|2.24p and the expressions (|2.26[) and (|2.28p . we verify that the 
expressions for the cross correlations obtained from formal substitution are correct. 

4.2 Imaging functionals using cross correlations 

We are now ready to present the imaging functionals, which consist in migrating the cross 
correlations of the recorded signals. The imaging functionals are designed according to the 
settings of receiver arrays. We consider two cases. 

Full aperture receiver array. The ideal case is when the receiver array spans the 
whole cross section, i.e., A = T>. Then given the data, for any pair of modes 4>j and (pi, we 
define 

C 3 \t):= [ I C(T,x 1 ,x 2 )(f> j (x 1 )(f>i(x2)dx 1 dx 2 . (4.8) 



I A J A 

Due to the orthoganality of {4>j}j=i ... n> the function l is the jl mode of the cross corre- 
lation function C. 

A search point for the reflector will be denoted as (a^i 3 ) where x s is its transversal 
coordinate and 2 s = L + z s is its axial coordinate. Equivalently, z s is the axial coordinate 
starting from the receiver array. We define the imaging functional Zfa as 

X FA (x s ,z s ) = l FA+ (x s ,z s )+l FA _(x s ,z s ), (4.9) 

N(w ) s 
l FA ±(x S ,z S ) = S2 PjPl(^o)^j{x S )(l>i(x S )C jl (±—(Pj+Pi)\. (4.10) 

3,1=1 

The choice of multiplication by /3j/3i is suggested from our analysis of the correlation func- 
tions in the next section. 

Limited aperture receiver array. A realistic situation is when the receiver array only 
covers part of the transversal section, i.e., A= [ai,a 2 ] for < a\ < a 2 < a. Consequently, 
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the exact jl component of the cross correlation cannot be extracted. In this case, we design 
the following imaging functional Zla as 

l LA (x s ,z s ) = l LA+ (x s ,z s )+l LA .(x s ,z s ), (4.11) 
1la±(x S ,z s ) = T W{ J 2 E M^)M^)M^(x S ) (4.12) 



"^o) 2 J a* 

:(A X1 + k(u ) 2 )(A X2 + k(u ) 2 )c( ± —(J3 q + fy) 



z s 



OJQ 

Again, these definitions of imaging functionals are based on the analysis of the correlation 
functions in the next section. We remark also that it is possible to show that when A = T), 
the second functional Xla is very close to 2fa and we have in fact: 

Xla±(x s ,z s )U =2? =T^ 7 ^2 E P]Pi(uo)M^)M^)C il (±—(Pi+Pi))- 
5 Resolution Analysis of the Imaging Functionals 

In this section, we analyze the imaging functionals proposed above to search the reflectors 
in the waveguide. Due to the random perturbations of the long section z G [0, L] of the 
waveguide, the values of the imaging functionals, which depend on the waveguide parameters 
through the data, are random. Hence, we analyze the mean of the imaging functional and 
show that it achieves its maximum at the reflector location. We also analyze how this mean 
decays from its maximum; this information provides the resolution of the proposed imaging 
functionals. 

We emphasize that the above observation on the mean of the imaging functional itself 
is not enough to claim that the functionals are effective, because it is not certain, a priori, 
that the one realization in practice is well reflected by the mean. Statistical stability (i.e., 
the analysis of the fluctuations of the imaging functionals) is needed to secure this claim. 
This will be investigated in the Section [6l 

5.1 The case of full aperture receiver array with time-harmonic sources 

We first consider the ideal case where the receiver array A spans the whole cross section, 
so the imaging functional Ip A is chosen. The key tool for analysis of the mean of Ifa is 
Proposition 13.11 

As shown by Proposition 13.1} when L is larger than the energy equipartition length 
L equ i p , the main contribution from terms of E{TJ^T^ re } comes from when j = m and I = n. 
Following this observation, we have the following limits for the cross correlations. 

E{C pp (t,x u x 2 )} -^-Yl W Mxi)Mx2)e-^ T , 

3=1 P] 

where 3>_x is defined by (|3.14p and we have also used the fact that converges to 1/N 
in the regime L S> L cqu i p . 
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Following the same lines, we have 

■ 2 <r r $-i(a g ) 
8 



(5.1) 

and 



E{C sp (r,X!,x 2 )} ^o^-ifo) £ -^^,(i 1 )^( a?r )^( a5r )^( a5l )e- i W«+ft)*e-** T . 

q,J = l J y 

(5.2) 

In fact, the migration imaging functional Xfa is designed from the above characterization 
of the cross correlation function. From the calculations before, we find 

where Sji is the Kronecker symbol. We also find 



E{Cg(r)} "" 0(J f" l( " s) ^^-(^)^(x r )e-^ + ^^-o-, 

Therefore, for a search point (x 8 ,? 5 ), with 5 s = L + z s , we have the following. From 
now on, A = 2ir/k(ujo) denotes the carrier wavelength. 

Proposition 5.1. If n> A, z s ,z r 3> A, and x s ,x T £ (0, a), i/ien 

E[Z FA (x s ,z s )] ~ ^lKe{[y' cos^ cose+ « sine )^] 2 }, (5.3) 

where we have introduced the normalized cross range offset £ = 27r(x r — x s )/A, and the 
normalized range offset fj = 2ir{z t — z s )/X. 

Therefore, the imaging functional Xfa works well to detect the point reflector (x r ,z r ). 
In particular we can see that both the range and cross-range resolutions are of the order of 
the wavelength A. 

Proof. Note that a>A means that N S> 1. Calculations show that 

/ \ -i N 



+ ^- 1 (x,) ^ (gS> ^ x ^ _ a,g^- l( x a ) ^ (a;S> zS . _^ )2) (5 4) 

where \& and <I>j are defined by f)3. 15|) and (|3.14p . The first term does not contain infor- 
mation about the reflector and can be viewed as background field. In fact, its contribution 
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is negligible because the fast oscillations in the complex exponential (z v ,z s S> A) aver- 
age out and there is an overall factor 1/N. For the second term, we can use an integral 
approximation for the sum in the continuum limit (N 3> 1). That is, 

,z ;x r ,z r ) = — 2^-sm(^— — )sm(^- — )e V a 
i=l 

N / r, / S\ ■ rv I S\ 

' J \7. zi—z- 
X 



1 s-^ f 2vr(x r -x b )j 2ir(x r +x b ) j \ i27T Ji-(±) 

j=i 

m-J [cos{ y)-cos( y)je l ^ ly a dy. 

Since the phase becomes zero when z s = z v , this function peaks at z s = z r and x s = x r . 
The integral considered above can be written as, with the second term neglected, 

_L f ^(vV^+iv) + ^(vV^'^dy = — C e^^^^dy 
2a J 2a J_ x 

7T 

= — I 2 cos ejfi ' cos6+i sin6) dO. 
2a J_il 

J 2 

Similarly, the sum in the third term can be approximated by 



*(x b , z b ; x r , -z r ) « - J \cos(^^ J -y) - cos (^ J X -y)) e V " ~dy- 

Note that this function does not have a peak comparable with the previous function. In fact, 
stationary phase calculation shows that it is of order 0(l/(ay k(z r + z s ))) where k = 2tt/X. 
The evaluation of E[Xfa— ] follows the same lines and we find a peak at z s = z r and 

X S = Xr- 

Finally we have 

N i 

which completes the proof of the proposition. □ 

Remark 5.2. To eliminate the (3j on the denominator, we multiplied by the exact f3j factor in 
constructing the imaging function. Alternatively, we can differentiate the cross correlation 
function to remove the denominator; see (I4.12p . If we do so, the above approximation will 
lead to the integral 

" cos 2 e «0«*»-«™< 



2 

in the square brackets of Eq. (|5.3p . This function should be compared with (|5.6p in the 
next section. 
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5.2 The case of limited aperture receiver array 



s 



Next, we consider the realistic setting where the receiver array only covers part of the 
transverse section, namely A = (01,02) and < a\ < a 2 < a. For a search point (x 
with 2 s = L + z s , we have the following result. 

Proposition 5.3. If 02 — a\ 3> A ; z s ,z r 3> a, andx s ,x r € (0, a), i/ien 



E[l LA (x s ,z s )] 



■Kuj^k o r (a2 — a\Y 
32aJ 



cos' 



Qgi(v cos sin 



(5.6) 



where we have introduced the normalized cross range offset £ 
normalized range offset fj = 2n{z r — z s )/X. 



2ir(x r — x s ) /X, and the 



This proposition gives the form of the point spread function of the mean imaging func- 
tional, that is the normalized form of the peak centered at the reflector location. The range 
and cross-range widths of the peak are the range and cross-range resolutions. Since the 
variables £ and fj are normalized with respect to the wavelength, we see that both the range 
and cross-range resolutions are of order of the wavelength, which is the diffraction limit. 
The form of the peak can can be seen in Figure [2] which plots the transverse and longitudinal 
shapes of the point spread function 



g(v) 



cos 2 9e^ sin6 d9 



7T 



(5.7) 



cos 



2 BS cosd d9 



cos 2 9 sm(fj cos 9)d9 ,(5. 



where we have used Formula 9.1.20 [T] to express h and g in terms of the Bessel function J\. 
h and g are even functions maximal at and stationary phase calculations also show that 



n(4) ~ 7T— 



sin(2g) 
T 3 



. , r?>i sin(2r/) 
V 



Proposition 15.31 also shows that the resolution of the mean imaging functional does not 
depend on the aperture of the array 02 — a\. This is a consequence of the waveguide 
geometry, since multiple reflections at the boundaries of the waveguide generate multiple 
replicas of the receiver array in the plane z = L, which gives an effective aperture that is 
large enough to reach the diffraction limit. 

Proof. Since the two functions above can be analyzed in the same manner, we focus on one 
of them, namely Xla+- As before, the cross correlation C can be treated component- wise. 
For the function Xla+, the main contribution comes from the primary-secondary component 
C ps , which we focus on for the moment. Let I+ ps denote this main term, i.e. 



1a 



ps 



% f N 

JA2 j',q'=l 



x(A Xl +k )(A X2 + k )C ps — (/3q> + Pj>),xi,x 2 )dxidx 2 . 
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Figure 2: Plots of the functions g(fj) and /i(£) which give the normalized form of the point 
spread function in the range direction (g(fj), left picture) and in the cross range direction 
right picture). 



Then from (|5,ip and the fact that A Xl (pj = —X 2 cj)j and k 2 = X 2 + /3? we verify that 



EfZ 4 



psj 



8N 2 



N 



Yl /3j/3q4>j(.Xr)4>j>(.X S )<t>q(x 



q,j,q',j'=i 

x / (f> j (xi)(f) j ,(x 1 )4> q (x2)4> q >(x2)dxidx2. 
J A 2 

If ^4. = (0,a), the last integral will be Sjj>5 qq ' and we are back in the case of full aperture 
array. In the current situation, this integral has to be dealt with more carefully. The above 
limit can be further written as 



v 



A 



(5.9) 



Introducing the difference index I = j — f , the double sum above can written in terms of f 
and I. The integral over x\ can be calculated explicitly as follows: 



(f>j'+i(xi)cj)ji(xi)dxi 



A 



a2 2 . (Z + Ovrxi . JW , 

- sm( ) sm( )ax\ 

a a a 



(12 



,irlx 



Q2 — a\ 



TT{2f + l)xi 

cosl ] — cos( jdxi 

a a 

f 7rl(a 2 + ai). . 7rl(a 2 - 



a 



cos 



2a 



-)sinc(- 



2a 



ir(2f + l)(a 2 + a 1 ) ir(2f + l)(a 2 -at) 
cos( )smc( ) 



2a 



2a 



For most j' , the first term above dominates. Hence, the number of indices Z's so that the 
above quantity is significant is roughly of order a/ (02 — a\). 
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Using the explicit expression (j)j(x) = \/2/a sin(jVx /a), we rewrite the function inside 
the brackets in (|5.9[) as 

f 1 ^Afij'+i . ir(j' + l)x T . ,7rj'x s . vr(j' + /)xi . ,-kj'xi 
J A N ^ a z a a a a 

e^A^-^e^'+^^dxL 

Using the explicit expression f3j = \Jk 2 — (jir/a) 2 and the fact that the total number of 
modes N is [a/(X/2)\ = [ak/n\ and assuming that ak/n is an integer for simplicity, we see 
that for / <C N, 

Using this expansion and some trigonometric identities, we approximate the integral above 
by 

-i ./ g . fc 2 jl l 

■ s^J X , ~ l ~ JV iV 2r 

sm( )e j 

a 

e^J^^dx!. 

Defined rj = z s — z r and £ = x s — x r . Substitution of some further trigonometric identities 
transform the above integral into 

/ "4 5>* { ( cos(^) cos(^) + S i„(^) sin(^)) co S (^) + R (j \(, l)\ 
J A a N 11 fl a a a J a J 

3 

(5.11) 

Here, I) consists of products of trigonometric functions, and for each product, one 

of the trigonometric functions is evaluated at 7rj'x v /a, irj'xi/a or irj'{x Y + x s )/a. Assuming 
that x r /a and x\/a are of order one, we observe that comparing with nj'^/a, which is 
explicitly written above and in which £ can be of order much smaller than one, irj'x r /a 
varies much faster as j' varies. In other words, Tij'x r /a can be viewed as a fast variable, 
and irj'^/a with ( < a is a slow variable. Similarly, in the expression of the complex 
potentials and fij'+i for fixed Z, the variable j'/N with N 1 involved is also a slow variable. 
Consequently, the contribution of R(j',£,l) to the sum over f is negligible comparing with 
those of the terms explicitly written. This two-scale analysis also shows that the above 
integral is much smaller for large £ comparing with £ <C a- 

We consider the regime TV 3> 1, and use the continuum approximation to rewrite the sum 
over j' as an integral with respect to the variable t = j'/N. We further assume 1 C I < JV, 
which is equivalent to say a ^> (02 — ai) 3> A; the second relation validates the linearization 
(|5.10p and allows us replacing f3j>+i by (3ji, while the first relation justifies the usage of 



A 



£ 

x 



4/3 



'j'+i 



a 2 N 



. ,TTJ X r 1llx t , ,7TJ X r . irlx r 

sm( ) cos( J + cos( J sin j 

a a a a 



2 .7Tj Hi, /kIx\ TTJXl /kIx\ 

sm ( Jcos( ) + cossm( sin ) 

a a a a 
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Poisson summation formula 

^Tcos^W^)-^^ = | E 

I meZ 



TT{X r + xx t 



irt z r 



+ 5C {X *- Xl) - -^=* + 2mn) + + - + 2mvr) 



a Vl - t 2 a 

7rt z r 



+ 2m7r) 



7T^2 a 



and 



E.7r/x r , .7rtxi N -t7r . - 
sm( )cos( )e v 1 -' 



2t 



' vr(x r + xi) vrt z r . 

7r^« +2m7r) 



-J( " (Xr " Xl) - -^=* + 2mn) - s C iXt " + —7=^— + 2mvr) 



a Vl - i 2 « 

7T(x r + ft Z, 

-0 ( h , h 2m7r) 



Note that we have used the fact that k/N = ir/a again. For each fixed m, let 5 a (x\,x T ; z T , m), 
a = 1, • • • ,4, denote the four Dirac distributions. Define also fj = 2tti]/X, £ = 2tt^/X; Then 
we have cos(nj'^/a) = cos(£i); the integral above becomes 



A 



^2 (Si + S 2 )(xi,x r ;z r ,m) 



+e i(VT=^fj+tO J- (5 3 + «5 4 )(xi, x T ; z T , m)dxi 



dt. 



If we extend the domain of t to (—1, 1), the above integral simplifies to 

kTT rl 



2a 2 



A 



^(5 3 + 8 4 )(x 1 ,x T ;z T ,m) 



dt. 



Integrate over xi first. The two Dirac distributions restrict the range of integration of t to 
the following intervals respectively: 



(ai + 2ma) — x T t 



(— ci2 + 2ma) — x r 



< 



< 



(a,2 + 2ma) — x T 



< 



(— ai + 2ma) — x T 



To interpret these conditions, imagine that the boundaries of the waveguide are two mirrors, 
then the array (01,02) has an image (—02, —ai) in the lower mirror; the two mirrors then 
generate a replica of images (ai + 2ma, 02 + 2ma) and (— 02 + 2ma, —01 + 2ma). If we define 
an angle 9 by 



9 = arctan 



t 



t€ (-1,1), 
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then the above restrictions of the Dirac measures can be restated as 

x T + a\ + 2ma x r + 02 + 2ma 
arctan < 9 < arctan , 

x r — 0.2 + 2ma x r — a\ + 2ma 
arctan < 9 < arctan . 

The first one restrict the angle 6 to those formed by the reflector and the array (a\, CI2) and 
the images of this array. The second set restrict the angle 9 to those formed by the reflector 
and the image array (—02, — a\) and the replicas. 

To analyze the resulting integral, we consider the simplest set-up where: (01,02) is 
centered in the cross section; the reflector is also centered in the cross range direction, i.e. 
x T = a/2. Further we assume the large distance regime z T S> a. 

In such a setting, with the notation a c = (01 + 02) /2 and w = (02 — a\)/2, the integral 
above boils down to 



1, /"arctan(ma+a c +to— x r )/z T 

— J2 cos 2 9e^ cos6+ ^ me U9. 

20 arctan(ma+a c — w— x T )/zi 



Since w <S z r , for each fixed m, the integral is over a very small angle section. Hence 
we can approximate the integral by the value at mean angle times the length of the angle 
section. The mean angle 6 m in the angle section is arctan ma/ z r . We further check that 



cos# m = l/y/l + (ma/z r ) 2 . Consequently, with a/z r set to Ax, the sum becomes 

k 2w z r y> exp{i(fj + mA^)(y / l + (jraAx) 2 ) -1 } ^ 
2a z r a ^ (1 + (mAx) 2 ) 2 X 



2a a (l + x 2 ) 2 2a 

Finally, let us verify that the primary-primary and secondary-primary cross correlations 
do not have significant contributions to the imaging functional Zla+- Let X +pp and I+ sp 
denote these two terms respectively. Similar to (|5.9|) . the expectation of X +sp converges, as 
e ->■ 0, to 



•? ^ 1 \ 1 N 



- £ f3 J 4> J (xrH f (x s )e-~ i ^ + ^' z > I ^(x^i^dxi ■ (5.12) 



2 



This function has the same form of (|5.9|) and can be analyzed in the same way. The key of 
these two functions is that the phase function in X +sp is a sum. As a result, the variable 77 in 
(|5.lip cannot be defined and have to be replaced by z v + z s which is of order a. This renders 
e %2-K^J\-{y /N) 2 (z T +z s )/x £ ag |. vai yi n g no matter how close z s is to z T . Due to the averaging of 
fast oscillations, there is no significant contribution from X +sp . 
For the primary-primary component, the analog to (|5.9p reads 



21 



Again, in the sum over f , only fast variables are involved. In the regime N 3> 1, the 
contribution of the function above is negligible. 

The term Xla- can be analyzed similarly Combining the main contributions in Ila+ 
and Xla-j we obtain the desired result. □ 



5.3 Imaging with Broadband Sources 

As we will see in Section El the imaging functionals are not statistical stable if the source 
is time-harmonic. Hence it is required to consider a broadband source (|4.6p . We show that 
the results obtained above for the means of imaging functionals apply to the broadband 
setting as well. 

Using Proposition l3.ll the main contribution of the two moment of mode coupling matrix 
at the same frequency comes from the case when j = m and / = n. Therefore, 

j=l J 

Similarly, for the primary-secondary field, we have 

iu^$_i(x s )cr r ^ 4>j{xi)<pj(x v )(/) q (x2)(l) q {x r ) 



Es«C ps (t, Xi , X2 ) j: 



e— >0 2tT J 

For the secondary-primary field, we have 

-iu%$-i(x s )a r ^ 4> q (xi)(j) q {x x )(j)j(x2)(l)j{x r ) 



E e a C ps (T, Xl ,x 2 ) — > u 2^ 



l ime -i(A+ft)^ e -^orJ_ /" |/ (M|V 
£-s>0 27T 7 



From these limits, we see that as long as z T , the distance between the reflector and the 
array is much smaller than e~ a , the integral in h above can be approximated by the energy 
of the source (square I? norm of /o). The rest parts of the limiting expectation of the 
cross correlation functions are exactly the same as the time-harmonic case. Consequently, 
the resolution analyses in the previous subsections based on the mean value of the cross- 
correlation migration imaging functionals remain the same. 



6 Stability Analysis of the Imaging Functionals 

The key tool is the following proposition which analyzes the asymptotic behavior of the 
fourth-order moment of the transmission coefficients in the limit e — > 0; see |10|. Section 
20.9.3] or HQ Section 8.4]. 
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Proposition 6.1. The expectation of four transmission coefficients at the same frequency 
has a limit as e — > 0. In the regime L 3> L cqu i p we have 

n(n+i) »/ 0'. = ( m > n ) = l ') = ( m '' n ') > 
jv ( ^ + i) *7 (i, = ("i, «) / (j'> = (m', n') , 

N{N+1) 0'> Z ) = ( m '' n ') ^ C/'j = ( m ' n ) > 

otherwise . 



Let a € (0, 2) and h 7^ 0. The expectation of four transmission coefficients at two frequencies 
ijj and uj + e a h has a limit as e — > 0. In the regime L 3> L oquip we have 

£>£e q uip J -L z/ (j,Z) = (m,n) and (j',/') = (m',n') , 



The previous sections shows that the mean of the imaging functional has a peak centered 
at the reflector location. The width of the peak is of the order of the wavelength. However, 
the imaging functional will give the reflector location only if it is statistically stable, that 
is to say, if the statistical fluctuations of the imaging functional are smaller than the mean 
amplitude of the peak. 

Time-harmonic case. We address the full aperture case in which the imaging func- 
tional is defined by (|4.10p . The mean of the imaging functional is (see (|5.4p ) 



E[X FA (x s , z s )} = ^d,_ l( x s ){K(^(x s , , s ; x r> z.f) +0{^)+ Q{j^-) }, (6.1) 



where &j and are defined by (|3.14p and (|3.15p . Here: 

- The term 0(1/ (kz s )) comes from contributions of the cross correlation of secondary (re- 
flected) and primary waves C ps and C sp that contain ^(x s , —z s ;x r ,z r ) or ^(x s , z s ;x r ,—z r ) 
(in such a case there are at least a product of two these terms, which gives the l/(kz s ) 
decay) . 

- The term 0(l/(NujQa r )) is the contribution of the cross correlation of primary waves C pp . 

The expression (|6.ip is valid provided N is large enough so that NoSqU^ 3> 1. Then it 
is true that the mean imaging functional is dominated by the first term in the right hand 
side, which is a peak centered at the reflector location. The mean amplitude of the peak at 
the reflector location is 

^pcak = ^-l(z s )$ (z r ) 2 . 



In the continuum approximation N 3> 1, we have (15. 5h and 



and therefore 



N 

t . 2 \— \ r,.27TX r j . JV»1 1 



Ppeak = (6.3) 
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The second moment of the imaging functional can be computed using Proposition 16.11 
in the regime E<1 and L S> L cqu i p . 

E[|X FA (x s ,, s )| 2 ] = ^L_\M[T FA (x s ,z S )}\ 2 

■ $ 2 i(xs) r $i(x S )$ _ i(;Er) |^ (a ,S^S ;;Cr ^ r) |2 



32(iV + 1 

+s( * (i ^s. IiiZi)1) + 0( _l )+0( i 



Here: 

- The term 0(l/(fcz s )) comes from contributions of the cross correlation of secondary (re- 
flected) and primary waves C ps and C sp that contain ^(a; s , — z s ;x r ,z r ) or ty(x s 

(in such a case there are at least a product of two these terms, which gives the l/(kz s ) 
decay) . 

- The term 0(1/ ' (Nuj 2 a T ) 2 ) is the contribution of the cross correlation of primary waves C pp 
that can be computes in a more quantitative way: 



°{ AT2 l A 2 ) = A r 2 4 4 2 frl(z S ) 2 + °( W*4 2 )- 

The variance of the imaging functional at the reflector location is 



-o( x r) 

In the continuum approximation iV ^ 1, we have (15. 5D . (16. 2p . and 

« io 



TV i 

^W = ^Eft™ 2 (— ^) - ;/ vT^dx = -, (6.4) 



and therefore 



Var(X FA (x r ,z r )) = Pp eak {- + - y}. 



To summarize: 

1) The typical amplitude of the fluctuations of the imaging functional for (x s , z s ) = (x r , z T ) 
(i.e. at the reflector location) is P pea k- 

2) The typical amplitude of the fluctuations of the imaging functional for |(x s ,z s ) — 
(x r ,z r )| 3> Ao (i.e. away from the reflector location) is -P pea k ((-WI^ S — ^l) 1 ^ 2 A (Ao|x s — 
Xl \) 3/2 + A§/(JVa r )). 

The second result shows that the relative fluctuations of the image due to the contri- 
butions of the primary cross correlation are of order X 2 ) /(Na I ). Therefore, provided the 
number of modes is large enough N 3> Ag/<7 r , these fluctuations are small. 

The first result shows that the amplitude of the peak at the reflector location has rel- 
ative fluctuations of order one. This is due to the fact that the reflector is illuminated by 
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a field whose amplitude is randomly spatially varying, so that the reflected energy is pro- 
portional to the squared amplitude of the primary field at the reflector location, which is 
a random quantity. This is the origin of the statistical instability in the time-harmonic case. 

Broadband case. We know that the frequency coherence radius in a waveguide with 
length L/e 2 is of the order of e 2 [101 Proposition 20.7] (or (TTJ Proposition 6.3]). As a 
result, as soon as a broadband source with a bandwidth larger than e 2 is used, then the 
field is the superposition of decorrelated frequency components. As a consequence the field 
is self- averaging in the time domain and we obtain that 

e 2 

Var(X FA (x r ,z r )) broadband ~ Var(X FA (x r , 2 r )) narrowband ^ 

where B 3> e 2 is the bandwidth of the source. We had seen that the use of broadband 
sources does not affect the resolution of the imaging functional but it ensures its statistical 
stability. Provided the bandwidth is larger than the frequency coherence radius, the typical 
amplitude of the fluctuations of the imaging functional is smaller than the amplitude -P pca k 
of the main peak at the reflector location, and therefore the reflector can be localized. 

7 Conclusions 

In this paper we have shown that the data recorded at a passive receiver array can be used 
in order to improve migration techniques in a random waveguide. The statistical stability of 
the imaging functional is ensured by the use of broadband sources. The resolution properties 
are ensured by the waveguide geometry: even when the receiver array does not span the 
whole cross section of the waveguide, the width of the point spread function of the imaging 
functional is of the order of the wavelength, provided the diameter of the array is larger 
than the wavelength. 

This paper has addressed the case of a two-dimensional waveguide with Dirichlet bound- 
ary conditions, but the conclusions should be qualitatively the same for fairly general situ- 
ations, when addressing three-dimensional waveguides, with Neumann, Dirichlet or mixed 
boundary conditions, with random fluctuations of the index of refraction or of the bound- 
aries as in [21 [T5] . 
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